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Abstract 

The quantum dynamics of coupled subsystems connected to a thermal bath 
is studied. In some of the earlier work the effect of intercenter coupling on the 
dissipative part was neglected. This is equivalent to a zeroth-order perturba- 
tive expansion of the damping term with respect to the intercenter coupling. 
It is shown numerically for two coupled harmonic oscillators that this treat- 
ment can lead to artifacts and a completely wrong description, for example, 
of a charge transfer processes even for very weak intercenter coupling. Here 
we perform a first-order treatment and show that these artifacts disappear. 
In addition, we demonstrate that the thermodynamic equilibrium population 
is almost reached even for strong intercenter coupling strength. 



PACS: 82.30.Fi, 82.20. Wt, 82.20.Xr 
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I. INTRODUCTION 



Quantum dynamics of complex molecules or molecules in a dissipative environment has 
attracted a lot of attention during the last years. One special kind of this problem is the 
electron transfer dynamics in or between molecules especially in solution |l]-f|. The bath- 
related relaxation can be described in a variety of ways. Among others these are the path 
integral methods [|||5| , the semi-group methods |||| , and the reduced density matrix (RDM) 
theory ||. The latter one has been especially successful in Redfield's formulation fT0|,[TT|l 
and is the topic of the present investigation. As usual, the master equation for the RDM is 
derived from the equation of the full system, i.e. relevant system plus bath, by tracing out 
the bath degrees of freedom. The main limitations of Redfield theory are the second-order 
perturbation treatment in system-bath coupling and the neglect of memory effects (Markov 
approximation). In addition Redfield suggested the use of the secular approximation. In 
this approximation it is assumed that every element of the RDM in eigenstate representation 
(ER) is coupled only to those elements that oscillate at the same frequency. In the present 
study we do not perform this additional approximation which could distort the correct time 
evolution in transfer problems |12|,p^]. 

To be rigorous in applying Redfield theory, the operators describing the time evolution 
have to be expressed in ER of the relevant system as has been done in the original papers 
limPI. F° r electron transfer this was performed in part of the literature (see for example 
|I1|-(P7|) while in another part of the literature [|TS|-|2"B|] diabatic (local) representations (DRs) 
have been used, which significantly reduces the numerical effort in many cases. In NMR 



literature [|n],|24|,^5j most people seem to use the ER while in quantum optics most people 
use DRs PE[p7|. Only recently the ER is used in quantum optics p8|-p0|. Here we focus 
on electron transfer systems, but the conclusions should also be applicable to problems in 
other areas. 

While in ER the damping term is evaluated exactly, in DR the influence of the coupling 
between the local subsystems on dissipation is neglected. As a consequence the relaxation 



terms do not lead to the proper thermal equilibrium of the coupled system [|,|3T],[32|. Only 
the thermal equilibrium of each separate subsystem is reached which can be quite different 
from the thermal equilibrium of the coupled system. It will be shown here that even for a 
very small intercenter coupling a completely wrong asymptotic value can be obtained. 

Although possibly leading to the wrong thermal equilibrium the local DR has advantages. 
For large problems it may be difficult to calculate the eigenstates of the system. These are 
not necessary in the DR. There one only needs the eigenstates of the subsystems. The 
quantum master equation can be implemented more efficiently in DR in many cases [|l8| ,|33l- 
35|. Moreover, almost all physical and chemical properties of transfer systems are expressed 
in the DR. For example, to determine the transfer rate one often calculates the population 
of the diabatic states and obtains the rate from their time evolution. To do so one has to 
switch back and forth between DR and ER all the time if the time evolution is determined 
in ER. 

Using the semi-group methodology and a simple model of two fermion sites, DR and 
ER have been compared already ||. We are interested in a more complicated system, i.e. 
a curve-crossing problem. The fact that we use a different relaxation mechanism should 
effect the findings only very little. Here we not only compare DR and ER but show how the 
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relaxation term in DR can be written more precisely for small intercenter coupling. 

The paper is organized as follows. The next section gives an introduction to the Redfield 
theory and, using the DR, presents a zeroth-order (DRO) and a first-order (DR1) pertur- 
bation expansion in the intercenter coupling. In the third section numerical examples are 
shown for two coupled harmonic oscillators. The DR and ER results are compared to each 
other and also to the improved local relaxation term derived here. The last section gives a 
short summary. Atomic units are used unless otherwise stated. 



II. INTERCENTER PERTURBATION EXPANSION WITHIN THE REDFIELD 

EQUATION 

In the RDM theory the full system is divided into a relevant system part and a heat 
bath. Therefore the total Hamiltonian consists of three terms - the system part H$, the 
bath part H B , and the system-bath interaction H^b'- 

H = H S + H B + H SB . (1) 

The RDM p is obtained from the density matrix of the full system by tracing out the degrees 
of freedom of the environment. This reduction together with a second-order perturbative 
treatment of H SB and the Markov approximation leads to the Redfield equation [P-|TT|, PE[ : 

p = -i[H s ,p]+np = £p. (2) 

In this equation 1Z denotes the Redfield tensor. If one assumes bilinear system-bath coupling 
with system part K and bath part $ 

#sb = K§ (3) 

one can take advantage of the following decomposition [06j |37||: 

p=-i[H s , P ] + {[Ap,K] + [K,pA% (4) 

Here K and A together hold the same information as the Redfield tensor TZ. The A operator 
can be written in the form 

oo 

A = J drmr)^(0))K l (-r) (5) 
o 

where K\—r) = e~ l Ht Ke is the operator K in the interaction representation. Assuming 
a quantum bath consisting of harmonic oscillators the time correlation function of the bath 
operator is given as |15| 



Of) = ($(r)$(0)) = J duJ(u)n(u)(e iwt + e^e~ iu}t ) . (6) 

o 

Here J(u) denotes the spectral density of the bath | 15| , n(uS) = (e^ w — the Bose-Einstein 
distribution, and (3 = l/(k B T) the inverse temperature. 
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The Hamiltonian Hs of the system we are interested in can be separated according to 



H S = H + V (7) 
where Hq is the sum of all uncoupled subsystem Hamiltonians Hq^ 

Ho = ^2 H o,n (8) 

n 

and V the coupling among them which is assumed to be small. Two canonical bases can be 
constructed for such a Hamiltonian. One consists of eigenf unctions of Hq. It is often called 
a local basis because these basis functions of the diabatic potential energy surfaces (PESs) 
are located at specific subsystems (centers). Latin indices such as \n) are used below to 
denote these DR basis states. The other basis diagonalizes the system Hamiltonian Ho,. So 
it consists of eigenstates of Hs and is called adiabatic basis. For these ER basis functions 
we use Greek indices such as \v). As discussed in the introduction Redfield theory is defined 
in ER but for transfer problems DRs have some conceptual and numerical advantages. 

Here we first calculate the dissipation in the DR for small intercenter coupling V . In this 
basis the matrix elements of A are given by 

oo oo 

(n\A\m) = J dujJ{uj)n{uj) J dr(e lulT + e^e~ lulT )(n\K\-T)\m) . (9) 



To evaluate the matrix element of K one has to use perturbation theory in V because the 
diabatic states \n) are not eigenstates of Hs but of Hq. Some details of the determination 
of (n\A\m) are given in the appendix. Using the expression for the correlation function in 
frequency space 

C(u) = 2tt[1 + n(w)] [J(oj) - J(-lu)] , (10) 

and denoting the transition frequency between diabatic states \m) and |n) by uj mn the final 
result can be written as 

(n\A\m) = -C(u mn )(n\K\m) 

1 (i\V\m\ 

^ j ^jm 

1 (n\V\i) 
- ^\K\m)^^[C{u mn )-C{u m ^] . (11) 

This first-order result DR1 can be split into a zeroth-order contribution DRO independent 
of V and a first-order contribution proportional to V. Taking the DRO term 

(n\A\m) = ^C(u mn )(n\K\m) (12) 

only is equivalent to a complete neglect of the influence of the intercenter coupling V on 
dissipation. This assumption has been used earlier |18|-|23|1 and is sometimes called the 



diabatic damping approximation ||38|| . In this approximation only the states |n) and |m) 
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contribute to the matrix element (n\A\m). In DR1 all states contribute to each of these 
matrix elements. As a consequence the spectral density of the bath is not only probed at 
the transitions of the uncoupled subsystems as in DRO but at many more frequencies. 

The ER result for the matrix elements of A can easily be deduced from the DR result by 
replacing the diabatic states by adiabatic ones and setting V = in Eq. (JTTJ) : 

(u\\\^) = ^C(u,J)(u\K\ij) . (13) 
This result is of course correct for arbitrary intercenter coupling strength. 



III. ELECTRON TRANSFER IN A TWO-CENTER SYSTEM 

In the following we direct our attention to electron transfer in an example system con- 
sisting of two charge localization centers considered to be excited electronic states. The 
PESs of the localization centers are assumed to be harmonic and are sketched in Fig. 1. For 
this example the Hamiltonian of the uncoupled system is given by 



U n + ( al a n + ^ ) io r 



(14) 



and the coupling by 



v = E E C 1 - S mn )v mn \mM)(nN\ . (15) 

m,n M.N 



The first index in each vector denotes the diabatic PES while the second labels the vibrational 
level. a n and at are the boson operators for the normal modes at center n and uo n are the 
eigenfrequencies of the oscillators. Bilinear system-bath coupling is assumed and the system 
part is given by the coordinate operator q 

if = g = EE (2u m M)- 1/2 (at, + a m ) \mM)(mN\ (16) 

m MN 

The mass of the system is denoted by A4. 

In the local DR the system part of the system-bath coupling reads 



(mM\K\nN) = {2u m My 1/2 5 mn (5 m +i,nVM + 1 + 5 M -i,nVm) . (17) 

In the DRO expansion (O) the system can emit or absorb only at intra-subsystem transition 
frequencies ujmn- The spectral density of the bath J (no) is effectively reduced to discrete 
values J(uS) = J2mlmS(uj — u) m ). The advantage of this approach is the scaling behavior of 
the CPU time with the number Af of basis functions which results from the simple structure 
of the A matrix (|I2D. As shown numerically |3~i|,p5] it scales like Af 2 ' 3 . This is far better 
than the Af 3 scaling of the DR1 approximation (10]). In DR1 the spectral density is probed 
at many more frequencies. One needs the full frequency dependence of J(u>) which we take 
to be of Ohmic form with exponential cut-off 

J(u) = rje^uje-"^. (18) 
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Here 9 denotes the step function and oo c the cut-off frequency. In this study all system 
oscillators have the same frequency uj\ (see Table |) and the cut-off frequency u> c is set equal 
to ui\. The normalization prefactor r\ is determined such that the spectral densities in DR 
and ER coincide at uj\. Eq. ( |18D together with Eq. (pi]) yields the full correlation function. 

If the system Hamiltonian H$ is diagonalized and the resulting ER basis is used to 
calculate the elements of the operators in Eq. (f|), there will be no longer any convenient 
structure in K or A, so that the full matrix-matrix multiplications are inevitable. For this 
reason the CPU time scales as A/" 3 , where M is the number of eigenstates of H§. There 
appears to be a minimal number A/o below which the diagonalization of Ho, fails or the 
completeness relation for \v) is violated. Nevertheless, the benefit of this choice is the exact 
treatment of the intercenter coupling. It is straightforward to obtain the matrices for p and 
K (see for example Ref. |jHH). 

An initial wave packet at center |n) is prepared by a 5-pulse excitation from the ground 
state |g) of the system 

PiMi7v(t = 0) = (lM|gO) ( g O|lA0 . (19) 

The pulse is chosen such that mainly the fourth and fifth vibrational level of the first (left) 
diabatic PES is populated. The motion of the initial wave packet along the coordinate q 
models the transfer between the centers. The parameters for our calculation are taken from 
the work of Kuhn et al. |2(J and are shown in Table |. Temperature is chosen as T = 295 K 



and the reduced mass of the system M. is set to 20 proton masses. The RDM is propagated 
in time and the occupation probabilities for each localization center are calculated by means 
of the partial trace: 

(20) 

M 

For the case of propagating in ER the RDM is transformed back to the DR in order to apply 
Eq. ©. 

In the following we compare the population dynamics in the two-center electron transfer 
system using three different intercenter coupling strengths V and four different configurations 
of the two harmonic PESs. The diabatic PESs and eigenenergies are shown in Fig. 1. 
Beginning our analysis with the weak coupling case v — v 12 = t>2i = O.lui it is expected 
that a perturbation expansion in V yields almost exact results. This is the reasoning why 
the DRO term, which is easy to implement, has been used in earlier work [I3-23|. 



In configuration (a) the eigenenergies of the two diabatic PESs are in resonance. For 
example, the first vibrational eigenenergy of the first center equals the third vibrational 
eigenenergy of the second center. It is important to note that in this configuration no vibra- 
tional level of the first center is below the crossing point of the two PESs. The calculations 
using ER and DRO as well as DR1 give almost identical results, see Fig. 2a. For long times 
DRO deviates a tiny bit. Redfield theory in ER is known to give the correct long-time limit 
(up to the Lamb shift). 

Configuration (b) differs from the first one by shifting the first PES up by cji/2. As 
shown in Fig. 2b the ER and DR1 results again agree perfectly. On the other hand, the 
DRO results are a little bit off already at early times and the equilibrium value departs from 
the correct value much more than in the first, on-resonance configuration. 
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Shifting the PESs further apart than in (a) yields configuration (c). The energy levels 
are again on-resonance but this time two vibrational levels of the first center are below 
the curve-crossing point, i.e. there is a barrier for low-energy parts of the wave packet. As 
shown in Fig. 2c DR1 and the ER results agree perfectly once more. The DRO results are 
terribly off. The long-time population of the first center which should vanish for the present 
configuration stays finite. 

If we increase the energy of the first PES by uj\/2 to obtain configuration (d) DRO fails 
again while DR1 gives correct results in comparison to the ER, see Fig. 2d. 

To understand the large difference between DRO and DR1 we have a closer look at the 
final result for the matrix elements of A, Eqs. ([11]) and (|T2|). The DRO contribution ( |i2|) is 



independent of the intercenter coupling V. The system part of the system-bath interaction 
K allows only for relaxation within each center. So there is no mechanism in the dissipative 
part which transfers population from one center to the other. This transfer has to be done by 
the coherent part of the master equation. But the coherent part cannot transfer components 
of the wave packet with energy below the crossing point of the PESs. As tunneling is mainly 
suppressed, those components of the wave packet cannot leave their center anymore although 
the corresponding PES might be quite high in energy. This results in the failure of DRO 
for the configurations with barrier: Parts of the wave packet get trapped in the two lowest 
levels of the left center. From Eq. (|TT|) one can explain why in the on-resonance case the 
DRO results are in better agreement with the correct results. In this configuration some of 
the DR1 terms are very small and so the DR1 correction is smaller. 

Now we discuss the medium coupling strength v = 0.5ui (see Fig. 3). The results for 
configurations (a) and (b), i.e. without barrier, look quite similar. In both cases the ER 
and DR1 results agree very well for short and long times. At intermediate times there is a 
small difference. The DRO results already deviate at short times and for long times there 
is too much population in the left (higher) center. For configurations (c) and (d), i.e. with 
barrier, again the ER and DR1 results coincide for small and long times. DRO is off already 
after rather short times and the long-time limit is again wrong. 

For the strong coupling v — u)\ (see Fig. 4) the behavior of the results is quite similar to 
the medium coupling. For configurations (a) and (b) the difference at intermediate times is 
a little larger, so is the deviation of the long-time DRO limit. For configurations (c) and (d) 
with barrier there is also a discrepancy for DR1 already at short times and the correct long- 
time limit is not reached exactly. But the disagreement is surprisingly small for the strong 
coupling. Overall DR1 still looks quite reasonable while the DRO results are completely off. 

IV. SUMMARY 

In addition to the approximations done in Redfield theory, i.e. second-order perturbation 
expansion in the system-bath coupling and Markov approximation, we have applied pertur- 
bation theory in the intercenter coupling. It has been shown for two coupled harmonic 
surfaces that the zeroth-order approximation DRO which is equivalent to the diabatic damp- 



ing approximation |38[ can yield wrong population dynamics even for very small intercenter 
coupling. These artifacts disappear using the first-order theory DR1. 

The scaling of DR1 is like J\f 3 not as A/" 2 ' 3 for DRO. This is of course a serious drawback 
of DR1. For configurations without barrier it seems to be possible to use DRO for weak to 
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medium intercenter coupling. This of course depends on the accuracy required especially 
for the long-time limit. In all other cases one should either use the exact ER or DR1. 
Although the first-order results are not exact for medium and strong intercenter coupling 
these calculations have at least two advantages. First of all, one does not need to calculate 
the eigenstates and energies of the full system Hamiltonian Hs- For small systems like two 
coupled harmonic surfaces using one reaction coordinate this calculation is of course easy. 
But if one wants to study larger systems like molecular wires JB|, |H| and/or multi-mode 



models |22|, |23], |33| this is no longer a trivial task. The second advantage is related to the 
fact that in all transfer problems one is mainly interested in properties which are defined in 
a local basis, e. g. the population in each subsystem in any moment in time. If one uses the 
ER one has always to transform back to the DR in order to calculate these properties. So 
for large-scale problems using a DR together with the first-order perturbation in V should 
be advantageous. 

In a sense the present study is an extension of the investigation performed by Davis 
et al. |J. They compared ER and DR for a two-site problem. Here we looked at a more 
general multilevel system and also calculated the first-order perturbation. In their model 
they do not have a reaction coordinate and therefore no barrier. Their findings correspond 
more to cases (a) and (b) in the previous section. Besides the agreement in the case of small 
intercenter coupling they also found good agreement in the high-temperature limit. Using 
our model this statement could not be confirmed for a general configuration, although there 
might be configurations where it is true. 

In Ref. |30f the authors followed a strategy different from the present work. They also 



studied two coupled harmonic oscillators modeling two coupled microcavities, but only one 
cavity was coupled to the thermal bath directly. This should not effect the questions studied 
here. With a transformation to uncoupled oscillators they effectively reduced the intercenter 
coupling to zero. The result |30| is then exact for arbitrary V. The disadvantage of this 
strategy is that it is not easy to extend to larger systems. The advantage of the presently 
developed first-order expansion in V is its general applicability to problems of any size. 
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APPENDIX: 

The purpose of this appendix is to show some more details for the evaluation of (n\A\m). 
To calculate 



(77 



K\-t)\m) = J2(n\e- im \i)(i\K\j)(j\e iHt \m) (Al) 



the operator identity 



-i(H 0+ v)t = e -iH t (l-ij dt'e u ' Ho Ve- lt ' (Ho+v) ^ , (A2) 



S 



which can easily be proven by multiplying both sides by e lHot and differentiating with respect 
to t, is used iteratively. It yields 



t 



(n\e- lHt \i) = (n\e- iHot [l - i J dt' e 11 ' Ho V e~ lt ' Ho ]\i) + 0{V 2 ) 

o 

t 

= e- iEit 5 ni - ie~ iEnt (n\V\i) J dt'e 1 ^-^' + 0(V 2 ) 

o 

= e~ iEit 6 ni - (e~ iEit - e~ iEnt ) + 0(V 2 ) (A3) 

assuming that E n ^ Ei. Here and in the following we only give the general expressions for the 
matrix elements. If a singularity can appear due to coinciding frequencies the appropriate 
expression can be obtained by taking the proper limit. 
Thus the matrix element (|A1|) is given by 

(n\K\-t)\m) = e iuJmnt (n\K\m) 

(j\V\m) 



- 52(i\K\m) {e iu3mit - e iuJmnt ) + 0(V 2 ) (A4) 

i ^ni 

This result is inserted into Eq. (Q). One has to evaluate integrals of the kind 

dte- et e- iUmnt = — (A5) 

^ — u nm — it 

which contain a convergence parameter e. Using the well known identity 

1 P 

lim = — =f ii5(x) (A6) 

e-*o x ± ie x 

one gets for the first term of the matrix element of A 

7T 

(n A \m) = w — [J(vmn) - J{-^mn)]{n\K\m) + (Lamb shift) + . . . 

The Lamb shift is the imaginary part of the matrix element of A and leads to an energy shift 
in the quantum master equation. This term is a small correction |4(],f41[] and is neglected in 
Redfield theory. The other terms of the matrix elements are calculated in the same fashion 
yielding 

7T 

(n\A\m) = - ■=- — [J{u mn ) - J(-u mn )](n\K\m) 

7T 



1 _ e -PWjr 



[J(uJ jn ) - J(-u jn )\ 
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(n|F|i) f 7T 
z|A|m)— < - [J(a; mra ) - J(-w mn )J 



7T 



[J(w mi ) - J(-w roi )]| (A7) 
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TABLES 



TABLE I. Parameters used for the ground state oscillator and the two excited state oscillators. 



Center |n) Configuration U n , eV Q n , A u> n , eV 

jg) O00 0.000 oX~ 

|1) 0.25 0.125 0.1 

1 2) a 0.05 0.238 0.1 

1 2) b 0.00 0.238 0.1 

|2) c 0.05 0.363 0.1 

|2) d 0.00 0.363 0.1 



13 



FIGURES 




FIG. 1. The four different configurations of the two diabatic harmonic potentials |1) and |2) as 
discussed in the text. Also included in the figures are the energy levels. 
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FIG. 2. Time evolution for small intercenter coupling and for the four different configurations. 

The results in ER are shown by the solid line while the results in diabatic basis are shown by dotted 

(zeroth-order) and dashed (first-order) lines. The results for ER and DR1 are indistinguishable for 

small intercenter coupling. Note the logarithmic time scale. 
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FIG. 3. Time evolution for medium intercenter coupling. 
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FIG. 4. Time evolution for strong intercenter coupling. 
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